/*
 * FieldInterpolation2D.cpp
 *
 *      Author: Jerome RD Soine and Christoph A Brand
 *      Institution: Schwarz goup, Institute for Theoretical Physics, Heidelberg University, Germany
 */


#include "../../../include/base/utility/FieldInterpolation2D.h"

FieldInterpolation2D::FieldInterpolation2D(std::vector<double> P1, std::vector<double> P2, unsigned int n_ref, std::string filename)
{
	Point<2> p1 (P1[0],P1[1]);
	Point<2> p2 (P2[0],P2[1]);

	triangulation = new Triangulation<2,2>();
	dof_handler = new DoFHandler<2,2> (*triangulation);
	fe = new FESystem<2> (FE_Q<2>(1), 2);

	generate_rectangular_mesh(p1,p2,n_ref);
	read_file(filename);
	interpolate_point_field(positions,displacements);
}


FieldInterpolation2D::FieldInterpolation2D(Point<2> p1, Point<2> p2, unsigned int n_ref, std::string filename)
{
	triangulation = new Triangulation<2,2>();
	dof_handler = new DoFHandler<2,2> (*triangulation);
	fe = new FESystem<2> (FE_Q<2>(1), 2);

	generate_rectangular_mesh(p1,p2,n_ref);
	read_file(filename);
	interpolate_point_field(positions,displacements);
}

FieldInterpolation2D::FieldInterpolation2D(Point<2> p1, Point<2> p2, unsigned int n_ref, SubstrateModel* substrate, vector<ForceTransmission*> fts)
{
	triangulation = new Triangulation<2,2>();
	dof_handler = new DoFHandler<2,2> (*triangulation);
	fe = new FESystem<2> (FE_Q<2>(1), 2);

	generate_rectangular_mesh(p1,p2,n_ref,fts);
	interpolate_fields(substrate, fts);
}


FieldInterpolation2D::~FieldInterpolation2D()
{
	delete data_container;
}

void FieldInterpolation2D::generate_vtk(string filename)
{
	output(filename);
}

void FieldInterpolation2D::read_file(string filename)
{
	data_container = new DataContainer(filename);

	unsigned int n_data_points = data_container->get_n_data_points();

	positions = new vector<Vertex*> (n_data_points);
	displacements.resize(n_data_points);

	for(unsigned int i=0;i<n_data_points;i++)
	{
		(*positions)[i] = new Vertex(data_container->get_position(i)[0],data_container->get_position(i)[1],0);
		Vertex tmp (data_container->get_displacement(i)[0],data_container->get_displacement(i)[1],0);
		displacements[i] = tmp;
	}
}

void FieldInterpolation2D::generate_rectangular_mesh(Point<2> p1, Point<2> p2,unsigned int n_refs)
{
	GridGenerator::hyper_rectangle(*triangulation,p1,p2,true);
	triangulation->refine_global(n_refs);
	dof_handler->distribute_dofs(*fe);
	solution.reinit(dof_handler->n_dofs());
}

void FieldInterpolation2D::generate_rectangular_mesh(Point<2> p1, Point<2> p2,unsigned int n_refs,vector<ForceTransmission*> fts)
{
	// Initialize global variables

		triangulation->set_mesh_smoothing(dealii::Triangulation<2,2>::smoothing_on_refinement);

		GridGenerator::hyper_rectangle(*triangulation,p1,p2,true);

		/*
		 * Use local refinement scheme.
		 */
		if(n_refs>=6)
		{
			triangulation->refine_global(6);
			n_refs=n_refs-6;
			for(unsigned int i=0;i<n_refs;i++){
					refine_based_on_traction_shapes(fts);
			}
		}
		else
		{
			triangulation->refine_global(n_refs);
		}
		dof_handler->distribute_dofs(*fe);
}


void FieldInterpolation2D::interpolate_point_field(std::vector<Vertex*>* positions, std::vector<Vertex> tractions)
{

	Vector<double> f (n_dofs());

	solution.reinit(n_dofs());


	DoFHandler<2>::active_cell_iterator cell;

		for(unsigned int i=0;i<positions->size();i++)
		{
			    std::cout << i ;
				Point<2> p((*positions)[i]->getX(), (*positions)[i]->getY());
				cout << "\t" << p[0] << "\t" << p[1] << "\n";
				cell = GridTools::find_active_cell_around_point<2>(*dof_handler,p);
				for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
				{
			      f(cell->vertex_dof_index(v,0))=tractions[i].getX();
			      f(cell->vertex_dof_index(v,1))=tractions[i].getY();
				}
		}

	solution = f;


}

void FieldInterpolation2D::interpolate_fields(SubstrateModel* substrate, vector<ForceTransmission*> fts)
{

	solution.reinit(n_dofs());
	force_solution.reinit(n_dofs());

	DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(),
										endc = dof_handler->end();

	for(;cell!=endc;++cell)
	{
		for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
		{
				Point<2> p = cell->vertex(v);
				vector<double> tmp_p (2,0);
				tmp_p[0]=p[0];
				tmp_p[1]=p[1];
				vector<double> vec = substrate->get_displacement_at(tmp_p);
			    solution(cell->vertex_dof_index(v,0))=vec[0];
			    solution(cell->vertex_dof_index(v,1))=vec[1];
			    for(unsigned int i=0;i<fts.size();i++)
			    {
			    	if(fts[i]->get_shape()->is_on_shape(tmp_p))
			    	{
			    		 double area = fts[i]->get_shape()->get_area();
			    		 force_solution(cell->vertex_dof_index(v,0)) = fts[i]->get_force()[0]/area;
			    		 force_solution(cell->vertex_dof_index(v,1)) = fts[i]->get_force()[1]/area;
			    	}
			    }
		}
	}

}

void FieldInterpolation2D::interpolate_force(vector<ForceTransmission*> fts)
{

	Vector<double> f (n_dofs());

	solution.reinit(n_dofs());

	DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(),
										endc = dof_handler->end();

	for(;cell!=endc;++cell)
	{
		for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
		{
			for(unsigned int i=0;i<fts.size();i++)
			{
				Point<2> p = cell->vertex(v);
				vector<double> tmp_p (2,0);
				tmp_p[0]=p[0];
				tmp_p[1]=p[1];
				if(fts[i]->get_shape()->is_on_shape(tmp_p))
				{
				vector<double> force = fts[i]->get_force();
			    f(cell->vertex_dof_index(v,0))=force[0];
			    f(cell->vertex_dof_index(v,1))=force[1];
				}
			}
		}
	}

	solution = f;
}

void FieldInterpolation2D::output(std::string filename)
{
  std::cout << "Writing solution vtk-file " << filename << "... \n";
  std::ofstream output (filename.c_str());
  DataOut<2> data_out;
  data_out.attach_dof_handler (*dof_handler);

  std::vector<std::string> solution_names;
  solution_names.push_back ("x_displacement");
  solution_names.push_back ("y_displacement");


  data_out.add_data_vector (solution, solution_names);
  data_out.build_patches ();
  data_out.write_vtk (output);
}

void FieldInterpolation2D::output_force(std::string filename)
{
  std::cout << "Writing solution vtk-file " << filename << "... \n";
  std::ofstream output (filename.c_str());
  DataOut<2> data_out;
  data_out.attach_dof_handler (*dof_handler);

  std::vector<std::string> solution_names;
  solution_names.push_back ("x_force");
  solution_names.push_back ("y_force");


  data_out.add_data_vector (force_solution, solution_names);
  data_out.build_patches ();
  data_out.write_vtk (output);
}

unsigned int FieldInterpolation2D::n_dofs()
{
  return dof_handler->n_dofs();
}

void FieldInterpolation2D::point_evaluation(Point<2> &p,Vector<double> &v)
{
	Functions::FEFieldFunction<2> fe_function(*dof_handler, solution);
	fe_function.vector_value(p,v);

}

void FieldInterpolation2D::refine_based_on_traction_shapes(vector<ForceTransmission*> force_trans)
{

	DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(),
												endc = dof_handler->end();

	for (; cell!=endc; ++cell)
	{
		 	 	   for(unsigned int i=0; i<force_trans.size();i++)
	    		   {
	    				vector<double> tmp = force_trans[i]->get_shape()->get_center();
	    				Point<2> p (tmp[0],tmp[1]);
	    				if(cell->point_inside(p))
	    					cell->set_refine_flag();
	    			}

	    			for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_face; ++v)
	    			{
	    				for(unsigned int i=0; i<force_trans.size();i++)
	    				{
	    					vector<double> vertex_point (2,0);
	    					vertex_point[0] = cell->vertex(v)[0];
	    					vertex_point[1] = cell->vertex(v)[1];
	    					if(force_trans[i]->get_shape()->is_on_shape(vertex_point))
	    					{
	    						cell->set_refine_flag();
	    					}
	    				}
	     			}

	}

	triangulation->execute_coarsening_and_refinement();

}
